library(HeadR)
rm(list=ls())
#
pref_reg <- "halton10k_2s"
#  
#Domestic (original C&G IV) ==== 
parnames <- fread("Output/Domestic/parameter_names.csv",header=FALSE)
parfiles <- dir("Output/Domestic/",pattern="^parameters_")
pars <- lapply(paste0("Output/Domestic/",parfiles),fread)
DL <- rbindlist(pars)
DW <- as.data.table(matrix(data=DL$V1,nrow=length(parnames$V1),ncol=length(parfiles)))
DW <- data.table(DW,var=parnames)
cols <- gsub(pattern="parameters_","",parfiles) |> gsub(pattern=".csv",replacement="") |> gsub(pattern="_domestic",replacement="")
setnames(DW,new=c(cols,"par_abbrev"))
setcolorder(DW,"par_abbrev")
#keeping only haltons
keepvars <- c("par_abbrev",pref_reg)
DW_domestic <- DW[, ..keepvars] # keep only preferred regression 
DW_domestic[, Domestic := ifelse(par_abbrev %contains% "s_" ,abs(get(..pref_reg)),get(..pref_reg))] # taking absolute value for sigmas
DW_domestic <- DW_domestic[, c("par_abbrev","Domestic")]
own_elas <- fread("Output/Domestic/own_elasticities_halton10k_domestic_2s.csv",header=FALSE) 
own_elas_domestic <- mean(own_elas$V1)
own_elas_domestic_1990 <- own_elas$V1[20]
rm(DW,DL,parfiles,parnames,pars,cols,own_elas)

#Domestic (Differentiation IV) ==== 
parnames <- fread("Output/Domestic/parameter_names.csv",header=FALSE)
parfiles <- dir("Output/DomesticDifIV/",pattern="^parameters_")
pars <- lapply(paste0("Output/DomesticDifIV/",parfiles),fread)
DL <- rbindlist(pars)
DW <- as.data.table(matrix(data=DL$V1,nrow=length(parnames$V1),ncol=length(parfiles)))
DW <- data.table(DW,var=parnames)
cols <- gsub(pattern="parameters_","",parfiles) |> gsub(pattern=".csv",replacement="") |> gsub(pattern="_domestic",replacement="")
setnames(DW,new=c(cols,"par_abbrev"))
setcolorder(DW,"par_abbrev")
#keeping only haltons (for dif iv we only have 1s)
pref_reg <- "halton10k_1s"
keepvars <- c("par_abbrev",pref_reg)
DW_domestic_difIV <- DW[, ..keepvars] # keep only preferred regression 
DW_domestic_difIV[, Domestic_difIV := ifelse(par_abbrev %contains% "s_" ,abs(get(..pref_reg)),get(..pref_reg))] # taking absolute value for sigmas
DW_domestic_difIV <- DW_domestic_difIV[, .(par_abbrev,Domestic_difIV)]
own_elas <- fread("Output/DomesticDifIV/own_elasticities_halton10k_domestic_1s.csv",header=FALSE) 
own_elas_domestic_difIV <- mean(own_elas$V1)
own_elas_domestic_difIV_1990 <- own_elas$V1[20]
rm(DW,DL,parfiles,parnames,pars,cols,own_elas)
pref_reg <- "halton10k_2s"

#Domestic NewIV ====
parnames <- fread("Output/Domestic/parameter_names.csv",header=FALSE) # same parameter_names
parfiles <- dir("Output/DomesticNewIV/",pattern="^parameters_")
pars <- lapply(paste0("Output/DomesticNewIV/",parfiles),fread)
DL <- rbindlist(pars)
DW <- as.data.table(matrix(data=DL$V1,nrow=length(parnames$V1),ncol=length(parfiles)))
DW <- data.table(DW,var=parnames)
cols <- gsub(pattern="parameters_","",parfiles) |> gsub(pattern=".csv",replacement="") |> gsub(pattern="_domestic",replacement="")
setnames(DW,new=c(cols,"par_abbrev"))
setcolorder(DW,"par_abbrev")
#keeping only haltons
keepvars <- c("par_abbrev",pref_reg)
DW_domestic_newIV <- DW[, ..keepvars] # keep only preferred regression 
DW_domestic_newIV[, Domestic_newIV := ifelse(par_abbrev %contains% "s_" ,abs(get(..pref_reg)),get(..pref_reg))] # taking absolute value for sigmas
DW_domestic_newIV <- DW_domestic_newIV[, c("par_abbrev","Domestic_newIV")]
own_elas <- fread("Output/DomesticNewIV/own_elasticities_halton10k_domestic_2s.csv",header=FALSE) 
own_elas_domestic_NewIV <- mean(own_elas$V1)
own_elas_domestic_NewIV_1990 <- own_elas$V1[20]
#standard errors for final regression
cov_diag <- fread("Output/DomesticNewIV/par_cov_halton10k_domestic_2s.csv",header=FALSE) 
cov_diag <- data.table(cov_diag,parnames)
setnames(cov_diag,new=c("cov_diagonal","par_abbrev"))
cov_diag[, se := sqrt(cov_diagonal/2217)] # reorder, merge, add column to table.
#
rm(DW,DL,parfiles,parnames,pars,cols,own_elas)


#Base regression ====
parnames <- fread("Output/Base/parameter_names.csv",header=FALSE)
parfiles <- dir("Output/Base/",pattern="^parameters_")
pars <- lapply(paste0("Output/Base/",parfiles),fread)
DL <- rbindlist(pars)
DW <- as.data.table(matrix(data=DL$V1,nrow=length(parnames$V1),ncol=length(parfiles)))
DW <- data.table(DW,var=parnames)
cols <- gsub(pattern="parameters_","",parfiles) |> gsub(pattern=".csv",replacement="") |> gsub(pattern="_domestic",replacement="")
setnames(DW,new=c(cols,"par_abbrev"))
setcolorder(DW,"par_abbrev")
DW[, Base := ifelse(par_abbrev %contains% "s_" ,abs(get(..pref_reg)),get(..pref_reg))] # taking absolute value for sigmas
DW_Base <- DW[, .(par_abbrev,Base)]
own_elas <- fread(paste0("Output/Base/own_elasticities_",pref_reg,".csv"),header=FALSE) 
own_elas_Base <- mean(own_elas$V1)
own_elas_Base_1990 <- own_elas$V1[20]
rm(DW,DL,parfiles,pars,cols,own_elas)


#NewIV  version of Base (no domestic) regression ====
parfiles <- dir("Output/NewIV/",pattern="^parameters_")
pars <- lapply(paste0("Output/NewIV/",parfiles),fread)
DL <- rbindlist(pars)
DW <- as.data.table(matrix(data=DL$V1,nrow=length(parnames$V1),ncol=length(parfiles)))
DW <- data.table(DW,var=parnames)
cols <- gsub(pattern="parameters_","",parfiles) |> gsub(pattern=".csv",replacement="") |> gsub(pattern="_domestic",replacement="")
setnames(DW,new=c(cols,"par_abbrev"))
setcolorder(DW,"par_abbrev")
DW[, NewIV := ifelse(par_abbrev %contains% "s_" ,abs(get(..pref_reg)),get(..pref_reg))] # taking absolute value for sigmas
DW_NewIV <- copy(DW)
DW_NewIV <- DW_NewIV[, c("par_abbrev","NewIV")]
own_elas <- fread(paste0("Output/NewIV/own_elasticities_",pref_reg,".csv"),header=FALSE) 
own_elas_NewIV <- mean(own_elas$V1)
own_elas_NewIV_1990 <- own_elas$V1[20]
rm(DW,DL,parfiles,parnames,pars,cols,own_elas)
#
#Merge all the results sets ====
DW <- merge_stata(DW_Base,DW_NewIV,by=c("par_abbrev"))
DW[, stata_merge := NULL]
DW <- merge_stata(DW,DW_domestic,by=c("par_abbrev"))
DW[, stata_merge := NULL]
DW <- merge_stata(DW,DW_domestic_difIV,by=c("par_abbrev"))
DW[, stata_merge := NULL]
DW <- merge_stata(DW,DW_domestic_newIV,by=c("par_abbrev"))
DW[, stata_merge := NULL]
DW <- merge_stata(DW,cov_diag,by=c("par_abbrev"))
DW[, stata_merge := NULL]
DW[, se_NewIV := paste0("(",round(se,3),")")]
#
fwrite(DW[,1:6],"Data/BLP99/pyblp_domestic_param.csv") # save csv of parameters for CF 
#
# Getting BLP99 parameters
param99 <- fread("Data/BLP99/published_param.csv") # BLP1995 parameters
setnames(param99,new=c("par_abbrev"), old=c("param_name"))
param99[, par_abbrev := sub(pattern = "const", replacement = "cons",par_abbrev)]
param99[, par_abbrev := sub(pattern = "sigma_", replacement = "s_",par_abbrev)]
param99[, par_abbrev := sub(pattern = "alpha_price", replacement = "b_price",par_abbrev)]
param99[, par_abbrev := sub(pattern = "demand_", replacement = "b_",par_abbrev)]
param99[, par_abbrev := sub(pattern = "supply_", replacement = "g_",par_abbrev)]
param99[par_abbrev=="b_price", pub_param := - pub_param]
DW <- merge_stata(param99,DW,by="par_abbrev")
DW[,c("pub_se","stata_merge") := NULL]
#
parnames <- fread("Output/Domestic/parameter_names_order.csv",header=TRUE)
DW <- merge_stata(DW,parnames,by=c("par_abbrev"))
DW[, stata_merge := NULL]
setorder(DW,order)
#
DW[,output:=texout(list(var,pub_param,Base,Domestic,Domestic_difIV),3)]
sink(file="Tables/CF_BLPdata/BLP_estimation.tex")
writeLines(DW$output)
sink()
#
own_elas <-data.table(var="mean own elas.",pub_param=-3.93,Base=own_elas_Base,NewIV=own_elas_NewIV,Domestic=own_elas_domestic,Domestic_difIV=own_elas_domestic_difIV,Domestic_NewIV=own_elas_domestic_NewIV,se_NewIV="")
own_elas_1990 <-data.table(var="mean own elas. 90",pub_param=-4.05,Base=own_elas_Base_1990,NewIV=own_elas_NewIV_1990,Domestic=own_elas_domestic_1990,Domestic_difIV=own_elas_domestic_difIV_1990,Domestic_NewIV=own_elas_domestic_NewIV_1990,se_NewIV="")
own_elas[,output:=texout(list(var,pub_param,Base,Domestic,Domestic_difIV),2)]
sink(file="Tables/CF_BLPdata/BLP_estimation_own_elas.tex")
writeLines(own_elas$output)
sink()

